A previous presentation showed that the Sundman transformation

$\frac{dt}{ds}=mr\sqrt{\frac{a^2(e^2-1)+2ar-r^2}{2m{r}^2(E-V)-L^2}}$

along with expressions for the energy

$E=\frac{1}{4e}[(1+e{)}^{2}V\left[a\right(1+e\left)\right]-(1-e{)}^{2}V\left[a\right(1-e\left)\right]]$

and the angular momentum

$L^2=\frac{m{a}^2(1-e^2{)}^2}{2e}[V\left[a\right(1+e\left)\right]-V\left[a\right(1-e\left)\right]]$

determines harmonic motion in the radial variable

$r=a(1-ecoss)$

for any spherically symmetric potential $V\left(r\right)$ in conjunction with the two quadratures

$t=\int dr\frac{mr}{\sqrt{2m{r}^2(E-V)-L^2}}$

$\phi =\int dr\frac{L}{r\sqrt{2m{r}^2(E-V)-L^2}}$

The motion is determined parametrically in an extension of the classical eccentric anomaly.

This presentation will show how to apply this solution in practice to power potentials of the form $V\left(r\right)=\alpha r^k$ . The point of adding the caveat “in practice” is that there can be a nonobvious difference between a mathematical object and its numerical evaluation.

First note that both energy and angular momentum squared are linear in the power potential, so that dimensioned quantities can be factored out:

$\begin{array}{c}E=\alpha a^k\frac{(1+e{)}^{k+2}-(1-e{)}^{k+2}}{4e}\\ L^2=m\alpha a^{k+2}(1-e^2{)}^2\frac{(1+e{)}^k-(1-e{)}^k}{2e}\end{array}$

Comparing these expressions with the two quadratures, it is clear the that second is independent of dimensioned quantities as it should be, being dimensionless itself. That means the three dimensioned quantities can be set equal to unity for its numerical evaluation. The first quadrature, on the other hand, has an overall dimensioned factor of $\sqrt{\frac{m}{\alpha a^{k-2}}}$ . Since the coupling constant has dimensions of energy over length to the power of the exponent in the potential, this works out to have the dimension of time as expected.

Taking $m=\alpha =a=1$ , a straightforward numerical evaluation of the second quadrature looks like this:

The calculation need only employ the parametrization by s in the limits of integration. Not only does this produce the same numerical result as rewriting the entire integral, it is also numerically much more stable. Further, moving the upper factor of angular momentum into the denominator of the integrand allows both positive and negative exponents to be handled together by canceling the negative coupling constant for the latter case.

The behavior of the integral is clearly a function of the periodic behavior of the radial variable, as can be seen by visualizing them together. This does not produce a result that can be used to generate an entire orbit, because the straightforward value never leaves the first Cartesian quadrant. What is needed is a constantly progressing value that can access all four quadrants.

A clue to how to proceed comes from thinking about the Kepler equation. The graph of that function is continuously ascending linear with sinusoidal variations about it. A similar behavior can be achieved here by flipping every other section of the function upwards modulus π, and shifting sections as appropriate for continuity. That produces this function:

This procedure can be justified in part by noting that the integral contains a square root, so that different portions of the full integral can be assigned different sheets of the square root function. Indefinite integrals also have a constant of integration which can be used to keep the function continuous.

The real justification is in being able to construct full rotating orbits that span all Cartesian quadrants. Since the semimajor axis is merely a scale on the overall orbits, keep constants as previously assigned. The parametric visualization of the orbits is then

This graphic takes a bit of time to process, because there is a significant amount of numerical math happening behind it.

Now consider explicit temporal evolution of the Cartesian coordinates. Again taking $m=\alpha =a=1$ for convenience, a straightforward numerical evaluation of the first quadrature looks like this:

The calculation again need only employ the parametrization by s in the limits of integration. Since there is no term proportional to coupling constant in the numerator, an absolute value is included in the denominator of the integrand to allow both positive and negative exponents to be handled together.

As before construct a continuous value by flipping every other section of the function upwards modulus π and shifting sections as appropriate for continuity:

This function is noticeably smoother than the angular quadrature. Now one can use both modified quadratures to show the individual Cartesian coordinates of the orbits. First the x-coordinate,

followed by the y-coordinate:

These last two graphics are both parametric plots. Interpolating the dependence of Cartesian variables on time can be done easily enough with a spline function if desired. Note that the quadrature for the temporal variable arises directly from the differential relation given at the outset, and as such is entirely equivalent to evaluating the latter.

As a final aside, consider the two quadratures from the start as instances of the integral

$\int dr\frac{r^n}{\sqrt{2m{r}^2(E-V)-L^2}}$

When the constants above are inserted in this general integral, it has an overall dimensioned factor of

$\frac{a^{n+1}}{\sqrt{m\alpha a^{k+2}}}=\frac{1}{\sqrt{m\alpha a^{k-2n}}}$

As before the coupling constant has dimensions of energy over length to the power of the exponent in the potential. In terms of dimensional analysis in units of mass, length and time, the overall factor becomes

$\frac{1}{\sqrt{M·\frac{M{L}^{2-k}}{T^2}·L^{k-2n}}}=\frac{T{L}^{n-1}}{M}$

For $n=1$ this is as expected time over mass, which explains the factor in the first quadrature. For $n=-1$ it is as expected the inverse units of angular momentum, again explaining the factor in the second quadrature. That opens up possibilities for other potentially interesting integral definitions of dynamic quantities.

Uploaded 2022.11.07 — Updated 2022.11.09 analyticphysics.com